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heat of formation of species i 
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Abstract 


This numerical study of the stability of laminar flames specifically 
addresses the dependence of stability on finite-rate chemistry with 
low activation energy and variable thermodynamic and transport 
properties. The calculations show that activation energy and details 
of chemistry play a minor role in altering the linear neutral stability 
results from asymptotic analysis. Variable specific heat makes a 
marginal change to the stability; variable transport properties, on the 
other hand, tend to substantially enhance the stability from a critical 
wave number of about 0,50 to 0,20, Also, the effects of variable 
properties tend to nullify the effects of nonunity Lewis number. When 
the Lewis number of a single species is different from unity, as is true 
in a hydrogen- air premixed flame, the stability results remain close 
to that of unity Lewis number. 

Introduction 

The linear stability of premixed flames has been studied by a large number of workers. 
Starting with Landau (1944) and Darrieus (1945), a number of workers, including Margolis, 
Matkowsky, Sivashinsky, and Clavin (see Clavin 1985 for a list of references), have contributed 
to this research. A scholarly review appears in Clavin (1985). Landau and Darrieus treated the 
one- dimensional flame as a thin discontinuity across which a density jump occurs. The dominant 
parameter that affects the stability is the disturbance wavelength in a direction normal to the 
direction of propagation. The investigations were conducted to find if such a lateral disturbance 
of wavelength w (or wave number k ) was amplified by the flow field, and the result was that 
the flame is unstable to all wave numbers. This result was puzzling for a long time because 
stable laminar flames could be set up in laboratories. Only recently have effects of nonnormal 
diffusion, curvature, and buoyancy on the flame stability been shown to reflect realistic features 
of stability. Studies by Sivashinsky, Matkowsky, and others included the effects of nonnormal 
diffusion (thermal and mass diffusion being unequal), first with low heat release and then with 
significant heat release. The analyses of Pelce and Clavin (1982), Matalon and Matkowsky 
(1982), and Frank el and Sivashinsky (1982) consider the limit of wave number tending to zero. 
The nondime nsionalizations used are different in these analyses, and the analysis of Pelce and 
Clavin is more general. One of the objectives of the three analyses conducted independently was 
to demonstrate the weak influence of viscosity in relation to conductivity and diffusivity. These 
studies obtain the dispersion relation (a result of the stability analysis giving the relationship 
between the growth rate and wave number) as to = a(q)k — b(q , Le)k , where to is the growth 
factor, k the wave number, and a and b are coefficients depending on the temperature ratio 
T^d/To and Le = Dpcp/n , the ratio of mass to thermal diffusivities. The expressions for a 
and b obtained by Matalon and Matkowsky (1982) as well as Frankel and Sivashinsky (1982) do 
not match (after appropriate transformations) even for a critical wave number. The numerical 
differences are not large in the range of parameters of interest, and the values predicted for 
a critical wave number are the same for Le = 1. Results based on systematic analysis and 
numerical integration of disturbance equations have been made by Jackson and Kapila (1984). 
Their numerical calculations have spanned the complete range of wave numbers, and they confirm 
the earlier results in the appropriate limits. They deduce from such analysis the influence 
of exothermicity and buoyancy on flame stability (Jackson and Kapila 1986). Increase in 
exothermicity is shown to destabilize the flame, and buoyancy is shown to stabilize the flame. 

All these studies are analytical in nature and have treated the high activation energy limit. 
In these studies the steady-state profiles for temperature, velocity, and mass fractions are 
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exponential in character and the reaction zone asymptotes to a plane. A typical plot from 
the calculations of Jackson and Kapila (1984) is shown in figure 1 for the largest parameter of 
exothermicity of Jackson and Kapila and without buoyancy. The parameter of exothermic ity 
corresponds to a flame temperature six times the cold reactant temperature. The abscissa in 

figure 1 is / = (jL— l)#. In this figure, m and 9 are proportional to the activation energy, and E 

for the finite-rate chemistry models is discussed later. The results obtained are for E /RT^ —> oo 
and for Lewis number not far from unity. It can be seen that at / = 0, corresponding to Le = 1, 
the unstable wave numbers are restricted to less than 0.36. And as Le increases, the range of 
unstable wave numbers increases. Values of Lewis numbers of some species go up to 2 in the case 
of H 2 ~air systems. For Le = 2.0 and E/RT^ — 4 for the H 2 ~air system, one obtains l — — 2 and 
*crit = 0-6. Thus the unstable wave number is increased substantially. At Le = 1, the unstable 
wavelengths are larger than 18 (= 27r/fc cr ^) times the flame thickness. The flame thickness is 
defined in these studies by 6^ = K r / Cp^ r p r u r , where K r is the reference conductivity, Cp^ r is the 
specific heat, p r is the reference density, and u r is the reference flame speed. For cold reference 
conditions, this gives 


0.052 W/m-K 

f ~ (1390 J /kg-K x 0.85 kg/ m 3 x 1.8 m/s) 

~ 0.0000244 m (0.0244 mm) 

For hot reference conditions, conductivity is almost 3.7 times the cold value and specific heat is 
25 percent higher. The product of density and flame speed is constant for the respective cold 
and hot values, and 6^ 0.072 mm. Therefore, the lateral wavelength causing instability is 18<5p 

about 1.3 mm. In the above results, the principal controlling factor is hydrodynamics. The role 
of diffusion seems significant only for Lewis number significantly different from unity. 

As stated earlier, in all the above analyses, the activation energy is treated as large. The 
overall activation energy has been estimated by Fenn and Cal cote (1957) to be 16 kcal/g-mole for 
the H 2 -air system and 28 to 30 kcal/g-mole for many stoichiometric hydrocarbon air systems. At 
typical flame temperatures of 2300 K, the activation parameter 9 — E / RT^ ^ 3.50 for H 2 -air 
and 6.1 to 6.5 for hydro car bon- air mixtures. Arguments concerning the validity of asymptotic 
analysis are made after estimating that 9 is between 10 to 20 (Clavin 1985). As noted previously, 
the perceived values of overall activation energy for equivalent single-step reaction are much 
smaller. Although this departure may still permit the validity of asymptotic analysis, there needs 
to be a demonstration of these aspects. Also, one may find that departures are small with regard 
to a few aspects, whereas for others, depending on the controlling phenomena, they are not. 

Present Work 

In the present work, the linear stability of flames is investigated numerically with particular 
reference to a stoichiometric H 2 ~air system by using a single-step finite reaction model. Two 
classes of reaction models, called “model A” and “model B,” are treated. Both are finite 
distributed reaction models. Model A is chosen because exact analytical solutions may be 
obtained of the steady state (or mean flow as it is called in stability analyses) and is primarily 
used to evaluate the effects of activation energy. Model B is in the line of classical single-step 
reaction models where numerical solutions are needed even for steady state. In this model the 
effects of variable thermodynamic and transport properties as well as the effects of diffusion are 
explored in detail. In the single-step reaction, there are four species, namely fuel (i — 1), oxidizer 
(i = 2), product ( i = 3), and inert ( i = 4). If the diffusion is modeled after the trace diffusion 
approximation (see Spalding 1957a and 1957b for details), one has four Lewis numbers, Le ? -, i — 1 
to 4, defined by Le ? - = pD^cp/hc, where Dj is the trace diffusion coefficient. Instead of letting Le ? - 
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vary through the held (something which is easily possible), different diffusion models are chosen 
to bring out the effects of diffusion explicitly. The first of these takes all the Lewis numbers to be 
equal to unity. This forms the reference case explored by the other investigators in asymptotic 
analysis. The second model takes all the Lewis numbers to be equal to 2.0. This value is chosen 
to represent the Lewis number of the hydrogen fuel in the environment of the other species. The 
third model takes Le^ = 2 and Le ? - = 1, i = 2, 3, 4. It is denoted by Le ? - = 2111. This is chosen 
because in a typical hydrogen-air flame, the Lewis numbers for species other than hydrogen are 
near unity. In addition, the sensitivity of the results to the accuracy of the steady-state pro hies 
is explored. 

Basic Equations 

The two-dimensional problem is set into an x-y Cartesian coordinate system, with the steady 
flame uniform in y and varying along x. A simple step reaction 

(Fuel + Oxidizer + Inert) — »■ (Product + Inert) 


is assumed. The conservation equations are as follows: 
Continuity: 

dp dpu dpv 

— + J —+ — = 0 
dt dx ^ dy 

^-momentum: 


p du du du dp f 4 d 2 u d 2 u 1 d 2 v 

~W +I,U sl +pv l3S 7+ ar r+ 3fcaJ 


y-momentum: 


p dv dv dv dp (4 d 2 v d 2 v 1 d 2 


dt 


+ pu H7 + pv + p \ + T 


dx dy dy \ 3 dx 2 dy 2 3 dx dy 


Energy: 


p dh s dh s dh s d ( dT\ d ( dTs ^ , , 

+ "“ ~ +pv W = Tz I s J + Ty I s ‘ 


dt 


dx 


°) Co m 

2 / l 


2=1 


Species conservation: 


p dYi , dYi , dYi d / dYA , 5 


dt 


dx 


dy dx \ dx J dy \ dy 


dYi 


+ P ull± + P v HI — — [ D- 1111 + — [D- HI l+l'" (* = 1,2,3) 

' r q __ 1 r 0 . . q __ \ 2 £> Q „ , 1 Q _ . I 2 £> Q . . , 1 1 V 5 5 / 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


pT _ p 
M ~ R 


( 6 ) 


where lu 1 ! 1 is the reaction rate of ith species. The primes denote the character of volumetric 
reaction rate. It is to be noted that for the single-step reaction, 


■ /// 
"l 


U) 


/// 


U) 


/// 


S+1 


(7) 


where s is the stoichiometric ratio and D 2 * represents the diffusion coefficients chosen in the 
present study to give the Lewis number of the ith species a desired value. The mass fraction of 
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the inert species (i = 4) is obtained from the identity that the four mass fractions must sum to 
unity. The equation of state assumes that pressure p is constant and is a good approximation 
for the stability study as well (Matalon and Matkowsky 1982). The momentum equations ignore 
the variation of viscosity in the held. This is assumed by noting the already demonstrated weak 
effects. (See Clavin 1985.) On the other hand, the variation of conductivity is accounted for 
because it is known to inhue nee the stability characteristics signihcantly. The equations are 
nondimensionalized as follows: 
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y 


tur 

= t 


— X 

k 

= y 


p_ 

Pr 

= P 

T 

T r 

= T 
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P_ 


Ur 

— V 

Kr 

— K 

pr 

= fi 


D iP 

(Dp) r 


Dip 


J 


( 8 ) 


The subscript r refers generally to cold upstream conditions, and 6^ is the hame thickness chosen 
as = K r / {pr u r c p ) r)' This implies that the Reynolds number based on the hame thickness is 


Re = 


p rUrSf 
Pr 


Kf* 

Pr c p,r 


1 

p7 


(9) 


Thus this choice of 6 ^ implies that the product of Reynolds number and Prandtl number of 
the hame is unity. Another nondimensional number which appears in the equations is Schmidt 
number 

§0 zz JL- 

(Dp)r 

which is set to unity to obtain (Dp) r . 

Dropping the bars over the symbols gives the following nondimensionalized equations: 

dp dpu dpv 


dt dx 


+ 


dy 


= 0 


p du du du dp ( 4 d^u d^u 1 d^v 

dt + P dx + P dy dx + 1 3 dx 2 dy 2 3 dx dy 


p dv 


dv 


di 


dt +PU d^ + PV Yy-~d^ + l 3 3? 


dp / 4 d^v d^v 1 d ^ 


+ 


+ 77 


dy ^ 3 dx dy 


Pr 


Pr 


p dh s dh s dh s d f dT\ d ( dT\ + h° Sf 

+ P U — + P v ^-Q^[ K —) + — t [ K 7^1 ~ 2^ 


dt 


dx 


dx 


at 


dy \ dy J ^ c p ^T r p r u r 


U)- 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 


P dYi 


OY; 


dY ;• d 


dYi 


+ pu — + pv — = — [ D iP — ) + — [ D iP — ) + 


d 


dy 


dt ' dx dy dx \ v dx 

pT — M s 

These are the basic equations we need to solve. 


dY; 


dy 


p f>U f 




(*' = 1,2,3) (14) 
(15) 
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Steady-State Equations 

Model A 

As discussed earlier, the one -dimensional model utilizes a formalism for which exact analytical 
solutions are available. The steady-state energy equation for model A is written as follows: 

d 2 r _ dr_ _ hctb^uj'" 
dx 2 dx p r u r (T ad - 1) 

where r — (T — ^/(T^ — 1) and h c ^ is the non dimensional heat of combustion. 

The thermodynamic and transport properties are taken as constant and the Lewis number is 
taken as unity. The species conservation equation is eliminated by adding to it the energy 
equation and obtaining the result of constant enthalpy through the flame. The continuity 
equation reduces to (pu) s = 1, and the momentum equation reduces to p — Constant after 
ignoring viscous and inertial terms. The choice of the reaction term is discussed in the section 
on solutions. 


Model B 

For model B the same approximations as for model A are utilized, and the energy equation 
is written similarly. The reaction equation and the rate expressions are 


2H 2 + 0 2 ^2H 2 0 (17) 

J" = A fP 3 Yl 2 Yo 2 exp (“Jr) - A bP Y E 2 0 ex P (jw) ( 18 ) 

The choice of the forward rate constants is discussed later. The backward rate constants are 
chosen to be consistent with the equilibrium constant for the reaction of equation (17). The 
resulting one-dimensional equations are solved by a code specifically developed for the purpose 
(Goyalet al. 1988). 

Stability Equations 

For stability analysis, the independent variables chosen are z — pu^ p, T, and Yp The 
various quantities are expanded around steady state, denoted by the subscript s, as 


p{x,y, t) = p s (x) + pf{x)<j>(y,t) " 
2 = z s (x) + Zf(x)<f>(y,t) 
u - u s (x) + Uf(x)<f>(y, t ) 

V — 0 + Vf(x)<j>(y ,t) > 

P = Ps(x) + pf(x)<f>(y,t) 

T = T s (x) +Tf(x)<j>(y,t) 

Y i = y s ,i( x ) + 


(19) 


In equation (19), z s (x) = pu s (x) = 1 because at steady state the equation of state gives constant 
mass flow through the flame. The disturbance function <j> is chosen as 

(f> = exp( — icut + iky) (i — \/ — l) (20) 


Note that iuj — l U r + and the sign of c or determines the stability of the flame. 
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The stability equations for x-momentum containing u f and u s must be expressed in terms of 
other quantities. Manipulations of the equation of state and the expression for z give 


Uf 


Ts_ , IL 

M s Z{ + M s 


U S = 


Ts_ 

M s 


and 


P f 


rj.Ms 

Tf W 


( 21 ) 


(22) 


where M s is the molecular weight variation at steady state. Perturbations of molecular weight 
are ignored in the analysis. 


Model A 

For model A, the following equations for zp pf , Vf, and Tf are obtained after the substi- 
tution of the expansions into equations (10) through (14) and the subsidiary relationships in 
equations (21) and (22) are used: 



(23) 


Pf - T 8 z'( + (t s ~ 2Tj) 


4 + 


'Pt s +Tg — T 7 ;) Zf + J(Ts)T[ = iLJZf 


a ! 1 2 -7 -wvf 

v { ~ v f ~ * 1 — 1 


u s 


(24) 

(25) 


T - T f ' + |j(T s ) - k 2 ] T { - -f Uf = 


. coTf 


(26) 


The primes denote derivatives with respect to x. This set of equations constitutes a seventh- 
order system. These are solved under conditions of zero values as x — »■ ±oo for all the 
variables zp pf, Vf, and Tf. Since the number of boundary conditions is eight, the problem 
is over determined and has nontrivial solutions for specific values of wave number k. Thus the 
problem becomes one of eigenvalue. 


Model B 


For model B, the following perturbation equations for z, x, p, and T are obtained when 
equations (19) and (20) are used in equations (10) through (14) and the subsidiary relationships 
in equations (21) are also used: 



— IUJ- 


T { M S 

Ti 


(27) 


Pf - Pr 




. M s 

f 

1 S 


(28) 


v i 


4Pr 


/ 

Vf 



kM’s 3 3 M s 

isf f - 


(29) 
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( 30 ) 


T " + Y + pwA r H ‘ 2 + i g Tt 


i ns 7 i io ns c 

1 hj + h; ^ 6r 


± u f T ^ C P rH .Ms^p r 

Tl^ ~r~Z^ ~^ J hJ Y j , f = ~ l ~ ^ 

k s “ Cn.rTr " /?ri/r ^ Ks TV k s 

%—\ 1 ’ j=l 


Y ff _ h^Y 4 - - I- 1 d{D^p) f f 1 ^f 

*-f *'W [ (D iP )Y (D iP ) s dT s s \ iY(D iP ) sPr u r 


J i,T J f + Y1 J i,j Y j, f 
3 = 1 


1 y, z 1 m?) , , M s 1 

D iPi Y '* ' tp iP ), il 1 ' 1 < T s (D iP ) s 't 


where 


c p — c p,s + T s 


j _ uu/ i 

^ “ dY; 


duj m 

_ uu/ l 

i,T~ Q T 


It is possible to relate the Jacobians J i j and J 2 * j- by invoking the stoichiometric relations 
between the reaction rates (eq. (7)) to obtain 


sJ l,j - J 2,j 


- (s + 1 )J ld - J SJ 


sJ IT — ^2 T 


- (s+ 1)J 1 T — J 3 T (34) 


In treating the variation of properties, it is assumed that all the dependence of the thermo- 
dynamic and transport properties in the flow field is described in terms of temperature alone. 
This is not entirely correct since there is some dependence of mass fractions of various species as 
well. But for premixed mixtures, it is reasonably accurate, certainly at a Lewis number of unity, 
where all the properties are described in terms of one progress variable, namely temperature. 
For nonunity Lewis numbers, the approximation implies that the extra dependence on mass 
fractions is ignored. The equation for pressure is in terms of Uf and its derivatives. They can be 
expressed in terms of Zf and Tf by using relationships in equations (21) to obtain an equation 
for pf given by 


/ Pr . // . o / / | A // | 4 rjiff nmJ 

Pf ~ — 4 u 8 Zf + 8n s z f + 4 u 8 z { + —T f - 8 T f 


AY) 2 


JM>) 2 

- i 772 T t + s ^zr T { 


9 9 TV i i / 1 / 74 / TV 

— Zk l u s Z£ — 3 k z — — + ikvr + 2uLzt + u s Zc H Tf hzMl — iujz f + iuj— (35) 

1 Ms f J ' S 1 f Ms f m] 1 T s V J 

The primes on various quantities represent derivatives with respect to x. The perturbation 
equations require u Sj derivatives u s (up to the second order), first derivatives of TV and Yp and 
Jacobians of reaction rate with respect to T and Y-. 

The order of the equations can be reduced further. One should normally solve three species 
conservation equations. However, the perturbation on the summation of mass fractions leads to 


y i,f + y 2,f + Y 3,f- 0 
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Hence it is sufficient to solve only two of the species conservation equations. The energy and 
species equations must be recast with the relationships given in equation (7) to get 


T{ + 


—c 


P 


+ 


2 dn 


k s k s dT 


Ts) 4 - T + 




K & pfVLy 


■JlT 7f 


1 h 

— fl c { 

Ks prUr 


(Jl,l - Jl, 3 )Y lf +(Ji , 2 ~ Jl, 3 )Y 2 f ~ ~T’z f = 

-I K s 


-,^T { (37) 

Is Ks 


Y 11 4 - 
y l,f + 


-k 2 + --T- Ylj 1 l -J 1 o) y, f - - 2 — Tig 2 ~J\ 3)^2 f 

(£V)s/V«r ’ ’ 6} \ ’ (Dip)s p r u r ’ 


1 


(Dip). 


-V 


i,f 


75 Y^ gZf + TT7 7 J 1 J'lf + 

D\p s ’ ( y D\p) s p r u r 


1 

(D iP ) s dT s 


Py! 


, 1 djpip) y! rpl _ 

+ D lPs dTs 1)5 f 


.Ms L_ 

(Di/9), 


-u>Y, 


l,f 


(38) 


Y n _|_ 

r 2,f + 


-*r 2 + — !_ ^(J 21 _ J 23 ) 

(Dip) s prur ’ ’ . 


y 2 s~ 


(D- 2 p) i 


%- 


(D‘2P) 


: Y 2,s z f + 


(D 2 P)s Pr u r 


(D 2 p)s prur 


J l,T T f + 


( J 2,l ~ J 2,3) Y 1,{ 

1 d{D 2 p) . , 

1 1 


( j ^2/ 9 )s d>T s 


s J l,f 


, 1 d P2^,/ M, 1 „ 

+ “ '^9 <JV — — ^ ” 77^ 


D<ips dT s 


T s (D2p)i 


(39) 


In obtaining these equations, the following equalities obtained from equations (34) have been 
used: 

J 2,2 “ J 2,3 = S ( J 1,2 ~ J l,3) J 2,l “ J 2,3 = S ( J 1,1 ~ J l,3) ( 40 ) 


^ct — 


h\ + + s (/^2 + ^2) — (1 + s) (^3 + ^3) 


Cp r T r 


(41) 


If the choice of Lewis numbers is such that Le 2 = Le 3 , then the two species equations can be 
related by sY^ f = — (s + 1)13 f . This equation can be used to reduce the number of species 
equations to one and to modify the energy equation. Such a modification gives 


dK 


T {' + -^ + 2k ^ T 's 4 + -r + -A 


l ct ' 


H 


c< 

v -+2k s 

k s dT s J L \ &s p r u r 

Cr) . M s Cp 


1 


J 


IT 1 i 


1 U S t T V C PrH - • 1V1 ST rri 

K S Ct p r U r 17 1)f Kg 1 T s /t s U f 


(42) 


Y [' { +{-k 2 + 1 


+ 


1 


Dips Pr u r 

h 


Ja Y lf + 


-1 1 dDip 

[{dW S + WT S dTs J ±Sll ’ f ~ WY) 


TlY'- 


hr Y u*t 


{&! p)s pr u r 


T ^ , 1 dD l Pvt rn! • M S 

d l,T^f + — — 77— 4| 0 T f — —no 


-Yi 


where 


(Dip) 8 dT s Ys f T s (D lP )s !> f 

J(T = JlA + s ^l,2 - (1 + s )^l,3 


(43) 

(44) 
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In this case the order of the equation to be solved is 9. If Le ? - = 1, one can eliminate the equation 
for Yi by combining it with the energy equation. The perturbation equation for enthalpy has 
zero for the solution. For this case, (Dip) s — and the temperature gradient terms associated 

c p 

with the transport properties are ignored. In this case the energy equation becomes 


rjiff P rri! 

J f - ~ 1 i 

Ks 


T f 


k 2 + — 




Ks prUr 


{Jcr ~ hct J\,t) 


Cp dT< 


S _ . Afg Cp 


Steady-State Solutions 

Model A 

Because equation (16) does not have the space coordinate explicitly, it is possible to reduce 

the order of the equation by defining q — ^ (Spalding 1957a and 1957b). One can then recast 
the equation as 

=“ a “" <46) 

where A represents the constants on the right-hand side of equation (46). Equation (46) has 
been analyzed in combustion literature. The reaction rate expression starts from exponentially 
small values near r = 0, peaks at some value of r depending on the activation energy, and 
goes to 0 at r = 1. Similarly, q is 0 both at r = 0 and 1 and is positive definite over the 
range r = 0 to 1 for the adiabatic case considered here. Based on these observations one can 
show that, for a class of pro hies q — r — r m , where m is a parameter, one obtains A = 1 
and uj nt — mr m (l — r m— ^). Reversing this argument, one can say that for this reaction rate 
expression (with m as a parameter), the solution for q is as stated earlier. One can integrate the 
equation for q and set out the steady-state solution as 


l 


T S 

— 1 + exp[— (m 

- 1 ) 

\x + c] >»- 


(47) 


T s 

— Ts ( -^ad 

- 1)+1 


(48) 



u s — 

T S 



(49) 

ii 

Eg ^3 

(Cad - 

»£= 

( T ad - 1 ) (r s 

-Tf) 

(50) 


■ m 

OJ - 

= ( 

1 - 

m— 1) 
T s ) 


(51) 

dJ" _ 

dT 

( r »i- 

l)m 2 T’” 

-1 

l-( 2 - 

_L ) T m - 1 

TO J 

(52) 


In these equations, the steady-state result that (pu) s = 1 along with the equation of state is 
used to obtain equation (49). In equation (47), c is chosen so that r = 0.5 at x — 0; this gives 

c = log(2 m_1 - 1) (53) 

The se solutions are coded and used in the solution of the stability equations. The choice of c 
has no effect except on the resolution of the eigen solutions. With the stability code using a 
grid distribution which allows a finer resolution at the center (x — 0) and increasingly coarse 
grid at x removed from this point, one expects better resolution by arranging the steady-state 
solution in this manner. The stability code utilizes its own grid distribution and computes the 
various quantities with the analytical expressions noted above. The parameter m characterizes 
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essentially the activation energy; that is, m = Constant(i?/i?T a( j). Typically values of m ^ 4 
to 6 imply high activation energy and m & 1.3 implies activation energies close to 16 kcal/mole 
of a range expected for a H 2 ~air system. This fact is based on the result that the reaction rate 
distribution (eq. (51)) has a peak at r s close to 0.5, a feature seen later even with a reaction 
rate distribution with r in the case of full chemistry. 

Model B 

The numerical solutions for the assumed reaction rate for the H 2 ~air system are obtained from 
an unsteady code developed for the purpose (Goyal et al. 1988). The code also generates the 
Jacobians of reaction rate with respect to temperature and mass fractions of species to be used in 
the stability code. These steady-state results are used in the stability calculations for nonunity 
Lewis number cases. The steady problem uses a uniform grid in z = f p dx. The grid is then 
transformed to the coordinate x — J dz / p. Then, the results of temperature, mass fractions, 
and Jacobians are interpolated with a cubic spline interpolation program into the grid required 
by the stability code. The temperature data then are spectrally differentiated by a Chebyshev 
polynomial fit to obtain the first and the second derivatives. These profiles were found to be 
jagged and nonsmooth. Consequently, it was decided to curve fit those data which needed to 
be differentiated. A Pade polynomial fit was used to describe the temperature distribution with 
x , and molecular weight, specific heat, conductivities, and diffusivities with temperature. These 
were then used in the stability code. 

Stability Solutions — Numerical Aspects 

The stability code used here was originally written for analyzing the stability of high-speed 
flows (Macaraeg, Streett, and Hussaini 1988). The perturbation equations are discretized by 
a spectral collocation technique using Chebyshev polynomials as basis functions. The code 
utilizes a staggered mesh to treat pressure. The resulting discretized equations are written in a 
generalized matrix eigenvalue problem and are solved with the standard library routine. (See 
Macaraeg, Streett, and Hussaini 1988.) 

Model A 

For model A, all the steady-state quantities were known in analytical form, and the 
calculations of the stability could be performed in a straightforward manner. The code utilized 
a grid stretching with the finest portion of the grids at x — 0. The region covered is from 
— oo to oo. It was therefore necessary to set a value for infinity. Several initial experiments 
suggested that infinity could be set at x — ±15. Sometimes, the eigenfunction could not be 
resolved accurately, since the decay was slow; for this reason, infinity was set at —20 to ±20. (It 
must be remembered that this x is already nondimensionalized by <5f.) Grid resolution studies 
were conducted and these showed that the results did not differ by more than 0.1 percent when 
the number of grid points exceeded 121. Most calculations utilized at least 121 grid points. An 
interesting aspect of the eigenfunction distribution was that pressure perturbations decayed the 
slowest toward the boundaries. Initial concerns regarding the effect on accuracy were resolved 
when it was determined that enhancing the boundaries and increasing the grid resolution did 
not affect the critical neutral wave number but altered the eigenfunctions marginally. 

Model B 

For model B, the range of infinity and the grid resolution used for model A were found to be 
valid. In the numerical results of steady flames, it was necessary to define a value of 6^. Although 
it would be possible to estimate the value from 6f = k r / p r u r Cpr , it was found convenient to 
assign a value to and with this value, obtain a consistent set of reference values. It should 
be remembered that the critical wave number, a result from the stability code, is actually a 
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nondimensional quantity, the nondimensionalizing parameter being One would expect the 
physical results obtained to be independent of the choice of This was ensured by varying the 
value of <5f and obtaining the critical wavelength for each case. 

Results and Discussion 

Model A 

Figure 2 shows the steady-state profile oftemperature, the first gradient, the second gradient, 
and the Jacobian. Most of the region of large change is restricted to a region — 6 < x < +6. 
The Jacobian varies significantly over the held, and the variation is different for the activation 
energy parameter m = 9 and 2. The variation is larger for m — 9 and smaller for m = 2. 
The calculations lead to a set of critical wave numbers (zero growth rate uq- = 0) for m, the 
activation energy parameter, as shown in table 1. The peaks of the eigenfunctions are shown in 
table 1(b) for the unperturbed (U) and the perturbed (P) cases, which are discussed later. It 
can be seen from the table that the critical wave number varies from 0.36 at high m to about 
0.40 at m ^ 1.3 corresponding to E 16 kcal/mole. This constitutes a 10-percent change which 
is not considered significant. The eigenfunctions are consistent with results from the asymptotic 
analysis. The imaginary part of zf, real part of vf , and imaginary parts of Tf and pf are zero. 
The other nonzero eigenfunctions are normalized by the peak of pf. 

In studies of stability with strong convection such as mixing or boundary layers, it was 
found that the mean profile exerts a significant influence on the stability characteristics. In 
order to determine the validity of this statement in the present context and to determine the 
features which affect stability significantly, subsidiary calculations were performed as follows. 
The initial profile of T f s {x) and T f s \x) was perturbed by a function Ax{xoq — x) / x^ sin ^r^f— 
chosen arbitrarily so that there would be fluctuations in the profile with zero at the boundaries 
x = 0 and x — Xqq. Figure 3 shows the plots of steady profiles. As can be seen, pro hies 
for both dT/dx and d^T/dx^ have considerable fluctuations. Table 1 shows and peak 

amplitudes of eigenfunctions for m — 9 and 2. The wave number k cv [ ^ is altered by no more 
than 3 percent, and the eigenfunctions are altered somewhat more but less than 10 percent. 
There are considerable fluctuations in the resulting eigenfunctions, largely those for pressure. 
These fluctuations do not seem to affect the overall result on stability. Thus the errors in 
temperature profile gradients seem to make little difference to the results of stability. The 
reason for this is that the instability is largely driven by hydrodynamics and details of the 
pro hie do not matter signihcantly. Figure 4 shows the eigenfunctions for both low and high 
activation energies, both perturbed and unperturbed cases. First, consider the unperturbed 
case. The structure of the eigenfunctions shows that their width is also from —6 to 6. It is only 
the eigenfunction for pressure that seems to decay slowly. For the lower activation energy, the 
temperature eigenfunction peak is larger than the pressure eigenfunction. This feature of the 
temperature eigenfunction having a peak higher than the pressure eigenfunction is seen in all 
the later calculations for model B. Between these unperturbed and perturbed cases the effect 
of disturbance is less severe for m — 9 than for m — 2. Calculations were made by changing 
the Jacobian by 5 percent from the nominal value. This results in a substantial change in the 
critical wave number of 20 percent. The features concerning the eigenfunctions look very similar 
and seem altered quantitatively to a small extent. Thus the stability is very sensitive to the 
Jacobiansbut quite insensitive to the details oftemperature profile gradients. 

The effects of Prandtl number have been discussed by earlier investigators (Clavin 1985) and 
were deduced to be insignificant . The results of the dependence of the critical wave number 
on Prandtl number are presented in table 2. The changes of k cv [ ^ near Pr = 1 are marginal. 
Only in the extreme case of Pr = 0.05 does the change of from that of Pr = 1 look 

substantial. A study that considered Pr — »■ 0 was conducted to determine if the viscous terms 
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could be neglected altogether. Two calculations were made by dropping the viscous terms in 
the momentum equations for u and v separately. Neglecting the viscous terms in the equation 
for u for m = 2 leads to a 10-percent reduction of k cv ^ from 0.391 to 0.356. Neglecting viscous 
terms in the equation for v does not lead to an acceptable solution satisfying the boundary 
conditions. This situation is inferred to be related to the neglect of the highest order derivatives 
in v (v n term), which is a typical singular perturbation problem. This is why the approach of 
obtaining the limiting solution of letting Pr — »■ 0 by retaining all the derivatives seems to lead 
to a physically consistent result. 

Model B 

Numerical calculations for a steady flame were performed for the stoichiometric H 2 ~air system 
with a single step reaction scheme, 2 H 2 + 02^ H 2 O, with frequency factor Aj of 1.1 x 10^, 
and an activation energy E of 16 kcal/mole. The choice of the parameters was based on the 
calculations of the stoichiometric flame structure with full chemistry (Bhashym et al. 1986). 
Figure 5 shows the plot of reaction rate of hydrogen with nondimensional temperature from 
such a calculation. The peak in the reaction rate occurs at T ^ 4.2, whereas the adiabatic 
temperature corresponds to T — 8.156. For Le = 1, the reaction rate expression becomes a 
function depending on temperature alone. Now, one can estimate E (or 9) from the plot of 
reaction rate with temperature. Such a calculation yields E & 16 kcal/mole. Such estimates are 
also available from earlier work (Fenn and Cal cote 1953). 

The steady flame speeds obtained from the steady-state calculations are 1.63 (Le 2 * = 1), 1.83 
(Le 2 * = 2), and 1.70 m/s (Le 2 * = 2111). The case Le^- = 2111 implies that the Lewis numbers 
for the four species 1, 2, 3, and 4 are 2, 1, 1, and 1. The results of the steady profiles and 
the eigenfunctions for the nominal case are shown in figure 6. The critical wave numbers for 
chosen values of 6^ are shown in table 3. As can be seen, the critical unstable wavelength is 
about 0.9 mm for the classical constant property case. For large activation energy, the critical 
wavelength would be about 1.05 mm (not shown in the tables). 

The calculations with variable properties show results which are interesting. Variable 
properties seem to act as a stabilizing influence, raising the unstable wavelength to as large 
as 1.88 mm. The property variation that has caused the change is deduced from the next two 
results. Variable specific heat alone seems to slightly destabilize the flame. But conductivity 
and diffusivity variation coupled through the Le^- = 1 assumption is the most stabilizing feature. 
It enhances the stability by a factor of 3. Clavin (1985) invoked the work of Clavin and Garcia 
(1983) and has indicated that the variable property effects can be taken into account by the use 
of thermal diffusivity at the hot condition rather than the unburnt condition. This effectively 
amounts to taking 6^ about 2 to 2.5 times higher than that estimated from the use of properties 
at unburnt condition. This effect then leads to enhanced stability. The results obtained in the 
current work are in conformity with the results of Clavin. The details can be understood by 
examining the results set out in figures 7 and 8. As can be seen from these figures, there are 
only weak differences in the profiles of the eigenfunctions, though the critical value of the wave 
number is significantly different between the constant and variable property case. 

Results of the kind described for model reaction were again established in the present case: 
(1) an increase of dT/dx by 1.5 changes the predicted critical wavelength by 2 percent, (2) a 
change in d^T/dx^ affects the results even less than a change in dT/dx , and (3) an increase 
of the Jacobian profile by 10 percent causes an increase in critical wave number of 25 percent 
(these results are not presented here). 

Once the range of infinity —15 to 15 and the scheme for interpolation were established, the 
approach to curve fit the steady-state quantities was abandoned in favor of numerically differ- 
entiating the temperature profile and using other interpolated quantities directly. Calculations 
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for constant properties turned out to be straightforward and gave the results within 1 percent 
of those from the curve fit noted above. The results for these calculations are presented in fig- 
ure 8. The second gradient alone seems to have fluctuations; but this does not seem to affect the 
eigenfunctions at all. The calculations for variable specific heat turned out to be more difficult 
to perform and needed better resolution. This was accomplished with 151 grid points. (The 
CPU times for the calculations of the eigen spectrum and the eigenfunctions on a CRAY-2 super 
computer were 71 s for 121 grid points and 120 s for 151 grid points for one case). 

Calculations have been made for two cases of nonunity Lewis number. In the first case, the 
Lewis number of all the species was 2. This corresponds to the conventional approach in which 
all Lewis numbers are equal. In the second case, the Lewis number of fuel alone is taken as 2.0 
and the Lewis numbers of the other species were set to unity. This follows from the calculations 
of H 2 -air either with full chemistry or single step chemistry with variable properties (see, for 
instance, Bhashyam et al. 1986) which show that Le ? - for H 2 is about 2 to 2.5, Le ? - for others 
is between 0.8 to 1.0. The steady-state profiles of Y 2 * s versus T and the Jacobi ans are shown 
in figure 9. The profile shapes for Le ? - = 2 show significant deviations from a linear profile. 
This is expected from simple analyses of the variation of nondimensional temperature with fuel 
mass fraction near a cold boundary (Spalding 1957a and 1957b). The profiles for Le ? - = 2111 
(Lewis number of various species in order of H 2 , O 2 , H 2 O, N 2 ), however, do not differ much 
from results with Le = 1. The differences in the results between Le ? - = 2111 and Le = 1 are 
caused by the diffusion terms. The stability results are summarized in table 4. The critical 
wavelength is typically 1.6 to 1.8 mm for the nonunity Lewis number cases. These values are 
only slightly smaller than for Le = 1. These observed features are a consequence of the fact that 
hydrodynamics controls stability and details of flame structure are less relevant to stability. 

Figure 10 shows the plots of the real and imaginary parts of the eigenvalues as a function of 
wave number k . It may be noted that l o r is less than 0 for the stable range shown in the figure. 
In all the cases except Le ? - = 2 (with variable properties), the imaginary part (uq-) is zero in 
the unstable range. The imaginary part being zero implies that the solution gets amplified in a 
nonoscillatory manner. The growth of the disturbance in time for any given unstable wavelength 
can be estimated from the results of figure 10. The time for doubling disturbance amplitude can 
be obtained from the disturbance equation (eq. (16)) as 


0.693 6 { 

uj f (k ) u f> 


(54) 


The time for doubling the amplitude scales like the characteristic time for the flame, with the 
coefficient typically being about 5 to 20. These values are relevant when making a full nonlinear 
simulation with a disturbance. 


Summary of Results 

The problem of the stability of laminar flames, particularly the H 2 -air system has been 
studied. The effects of finite-rate kinetics and variable thermodynamic and transport properties 
are explored. The perturbation equations are spectrally discretized and numerically solved 
to obtain the eigenvalues and the corresponding eigenfunctions. These calculations show the 
following results: 

1. The effect of finite activation energy on the critical wavelength is not significant. Reduction 
of the activation energy to values corresponding to the H 2 -air system reduces the critical 
wavelength by about 10 percent. 

2. Variable transport properties enhance the stability and enhance the critical wavelength by 
a factor of 2 to 2.5. 
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3. Results for realistic parameters show that the critical unstable wavelength for a stoichio- 
metric H 2 ~air mixture is about 1.6-1. 8 mm. 


NASA Langley Research Center 

Hampton, VA 23665-5225 

November 19, 1991 

References 

Bhashyam, A. T.; Deshpande, S. M.; Mukunda, H. S.; and Goyal, G. 1986: A Novel Operator-Spitting Technique 
for One -Dimensional Laminar Premixed Flames. Combust . Sci. & TechnoL, vol. 46, nos. 3-6, pp. 223-248. 

Clavin, Paul 1985: Dynamic Behavior of Premixed Flame Fronts in Laminar and Turbulent Flows. Prog . Energy & 
Combust . Sci., vol. 11, no. 1, pp. 1-59. 

Clavin, P.; and Garcia, P. 1983: The Influence of the Temperature Dependence of Diffusivities on the Dynamics of 
Flame Fronts. J. Mec. Theor . & AppL, vol. 2, no. 2, pp. 245-263. 

Darrieus, G. 1945: Propagation d’un Front de Flamme. Paper presented at Congres de Mecanique Appliquee (Paris). 

Fenn, John B.; and Calcote, Hartwell F. 1953: Activation Energies in High Temperature Combustion. Fourth 
Symposium ( International ) on Combustion ( Combustion and Detonation Waves), Williams A Wilkins Co., 
pp. 231-239. 

Franket M. L.; and Sivashinsky, G. I. 1982: The Effect of Viscosity on Hydrodynamic Stability of a Plane Flame 
Front. Combust . Sci. & TechnoL , vol. 29, nos. 3-6, pp. 2 07-224. 

Goyal, G.; Paul, P. J.; Mukunda, H. S.; and Deshpande, S. M. 1988: Time Dependent O per at or- Split and Unsplit 

Schemes for One Dimensional Pre mixed Flames. Combust . Sci. & TechnoL , vol. 60, pp. 167-189. 

Jackson, T. L.; and Kapila, A. K. 1984: Effect of Thermal Expansion on the Stability of Plane, Freely Propagating 
Flames. Combust . Sci. & TechnoL , vol 41, nos. 3-4, pp. 191-201. 

Jackson, T. L.; and Kapila, A. K. 1986: Effect of Thermal Expansion on the Stability of Plane, Freely Propagating 
Flames. Part II: Incorporation of Gravity and Heat Loss. Combust . Sci. & TechnoL , vol. 49, pp. 305-317. 

Landau, L. 1944: On the Theory of Sbw Combustion. Acta Physicochim ., vol. 19, p. 77. 

Macaraeg, Michele G.; Streett, Craig L.; and Hussaini, M. Y. 1988: A Spectral Collocation Solution to the 

Compressible Stability Eigenvalue Problem. NASA TP-2858. 

Matalon, M.; and Matkowsky, B. J. 1982: Flames as Gasdynamic Discontinuities. J. Fluid Mech., vol. 124, 

pp. 239-259. 

Pelce, Pierre; and Clavin, Paul 1982: Influence of Hydrodynamics and Diffusion Upon the Stability Limits of Laminar 
Premixed Flames. J. Fluid Mech., vol. 124, pp. 219-237. 

Spalding, D. B. 1957a: I. Predicting the Laminar Flame Speed in Gases With Temperature -Explicit Reaction Rates. 
Combust . & Flame , vol. I, pp. 287-295. 

Spalding, D. B. 1957b: II. One -Dimensional Laminar Flame Theory for Temperature- Explicit Reaction Rates. 
Combust . & Flame , vol. I, pp. 296-307. 


14 



Table 1. Model A 


(a) Critical wave number k cv - lt (b) Peak amplitudes of eigenfunctions 


m 

Case 

Re(V) 

Im(T) 

Re(T) 

Re(p) 

9 

Unperturbed 

2.9 

17.5 

74.0 

210.0 

9 

Perturbed 

3.3 

19.2 

79.1 

225.6 

2 

Unperturbed 

8.2 

28.1 

78.3 

58.3 

2 

Perturbed 

9.1 

29.3 

79.7 

62.5 


m 

&crit for — 

Unperturbed 

case 

Perturbed 

case 

9 

0.362 

0.3704 

4 

0.375 

0.382 

2 

0.391 

0.407 

1.3 

0.401 

0.411 


Table 2. Model B 


Pr 

k crit for 

m of — 

9 

2 

1.0 

0.391 

0.362 

0.7 

0.377 

0.374 

0.1 

0.305 

0.331 

0.05 

0.300 

0.301 


Table 3. Model B for Le ? ; = 1 and 9 — 3.5 

(a) Critical wave number & cr p and (b) Peak amplitudes of eigenfunctions 


Case 

h • 

^ crit 

mm 


Constant properties 

0.42 

0.06 

0.9 

Variable properties 

0.20 

0.06 

1.88 

Variable c p 

0.43 

0.06 

0.88 

Variable tz and Dp 

0.18 

0.06 

2.09 


Case 

Re(V) 

Im(V) 

Re(T) 

R e(p) 

Constant properties 

BB 

26.0 

83.8 

64.5 

Variable properties 

HI 

25.0 

74.4 

66.5 


Table 4. Model B for Le ?: = 2 and 2111 


(a) Critical wave number & cr p and 


(b) Peak amplitudes of eigenfunctions 


Case 

Le; 

^crit 

St, 

mm 

Him 

Constant properties 

2 

0.50 

0.06 

0.75 

Variable properties 

2 

0.24 

0.06 

1.57 

Constant properties 

2111 

0.48 

0.06 

0.785 

Variable properties 

2111 

0.21 

0.06 

1.79 


Case 

MB 

Re(«) 

Im (r>) 

Re(T) 

Re(p) 

Constant properties 

2 

9.5 

23.0 

86.4 

60.4 

Variable properties 

2 

6.2 

21.0 

84.4 

64.5 
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